seem (version 1.0)

Model eval and sensitivity: Model evaluation and sensitivity analysis

Description

Function model.eval provides calculation of errors between observed or reference values and simulated values.

Function sens provides sensitivity analysis on one parameter using six metrics. Function msens conducts sensitivity analysis on combinations of several parameters. These are combined by full-factorial (using mway.fact) or stratified sampling by Latin Hyper-cubes (using mwy.samp).

Usage

model.eval(t, Xobs, Xmod, fileout, pdfout=F)
sens(t.X, param, fileout, pdfout=F)
msens(t.X, param, fileout, pdfout=F, resp=F)
mway.fact(pnom,plab,perc.v,levels)
mway.samp(pnom,plab,perc.v,runs)

Arguments

t

times to evaluate output

Xobs

observed values

Xmod

simulated values

fileout

name of output file

pdfout

logical to decide to produce PDF output file

t.X

output of multiple runs for sensitivity

param

list of parameter variation for sensitivity, nominal, label, values, and a logical variable fact for factorial or Monte Carlo analysis

resp

logical variable to decide upon producing values of response function

pnom

array of nominal parameter values

plab

array of parameter labels

perc.v

percent variation around nominal

levels

levels for full factorial analysis

runs

number of runs for Monte Carlo

Value

• For model.eval
• result
  Table of individual errors
• label.result.tot.avg
  Labels to identify the columns of total error
• result.tot.avg
  Total errors
• For sens
• val.val
  metric values corresponding to param values
• perc.perc
  percent change in the metric to percent change in the parameter value
• perc.ratio
  percent ratio change in the metric to percent change in the parameter
• reg.slope.R2
  slope and R-square of regression lines for sensitivity analysis of Monte Carlo runs
• For mway.fact
• plab
  array of parameter labels
• pval
  combined parameter values
• pnom
  array of nominal values
• fact
  echo of fact argument
• pv
  table of parameter values
• For mway.samp
• plab
  array of parameter labels
• pval
  combined parameter values
• pnom
  array of nominal values
• fact
  echo of fact argument
• For msens
• Ma
  values of parameter values and six metrics for each run
• Mp
  percent change of parameter and six metrics for each run
• x
  When resp=T, values to make these response maps. This list has six elements, one for each metric. The first four of each set of six correspond to the metrics. The other two refer to the difference with respect to nominal.
• AOV.pvalue
  when fact=T, p-values of the ANOVA to examine whether there is a significant difference in the response of the metrics
• Friedman.pvalue
  when fact=T, p-values of the nonparametric test to examine whether there is a significant difference in the response of the metrics
• Eff
  when fact=T, average effect for each parameter on all metrics
• Eff.Int
  when fact=T, effect for each parameter in interaction with the others on all metrics
• reg.slope.R2
  When fact=F, slope and R-square of regression lines for sensitivity analysis of Monte Carlo runs

Details

Function model.eval calculates absolute, relative and square errors or each pair of observed and simulated data. It also provides total absolute, relative, and RMS.

Function sens operates on output from multiple runs made by varying a parameter and examining the change of one or several metrics of the output. These all apply to X(t) for the entire simulation run; that is, a metric summarizes the entire run in just one value. Therefore, we can analyze the effect of changing the parameter on the value of the metric and plot x-y graphs for visualization.

The six metrics are: g1() is the maximum value of Xi for run i; g2() is the minimum value of Xi for run i; g3() is the value of Xi for run i; g4() is the average of the difference Xi-Xn for run i; g5() is the RMS of the difference Xi-Xn for run i; and g6() is the maximum of the absolute value of the difference Xi-Xn for run i. Note that the last two metrics always yield positive values. Their sign is inverted for values corresponding to the parameter values below the nominal.

For function sens: often, we perform an odd number of runs where the nominal value corresponds to the center value in the set of parameter values, and we vary by a percent. However, the nominal can be anywhere as dictated by the component pnom in the param argument.

Function msens: The pairwise combinations for more than two parameters employ the function factorial and the function tapply. The latter applies the mean to the response using the parameters as factors and produces tables x and xp, which are the basis to calculate the effects and maps of the response. The average effects are computed from the tables xp. Multiple regression is performed using lm for each metric versus the relative parameter values, ps, which had been previously standardized so that we can compare coefficients across the parameters. We extract the coefficients and the R2 from the summary of lm, which is a list. Then, we use the coefficients to calculate the predicted metrics in order to employ these later to graph the scatter plots of metric calculated versus metric predicted. Finally, all information is packed in a data.frame to be returned by the function. Interaction plots are produced using the function interaction.plots. Response maps are generated by means of the function image and the overlapping isolines are generated using the function contour with the option add=T.

The function mway-fact helps determine a full factorial design. For example, we can design the factorial of three levels at 50 percent. Function mway.fact uses two R functions, gl and factor. The output of mway.fact is used as an argument to simulation functions (e.g, sim) and to msens.

Function mway.samp is for stratified sampling of more than one parameter that varies at a time. It performs Latin hypercubes to obtain n numbers by sampling using a uniform distribution from a min to a max value determined from the desired percent variation. It is limited to the uniform distribution but can be extended to sample other distributions. The output of mway.samp is used as an argument to simulation functions (e.g, sim) and to msens.

References

Acevedo M.F. 2012. Simulation of Ecological and Environmental Models. CRC Press.

Kirchner, T.B. 1992. TIME-ZERO The Integrated Modeling Environment. Reference Manual. Quaternary Software.

See Also

Simulation functions sim.comp, sim.mruns

Examples

# model evaluation
# read observed data from file
t.Xobs <- matrix(scan("chp5/expodecay-data.txt", skip=1), ncol=2, byrow=TRUE)
# simulate
model <- list(f=expon)
# assign values for calculating errors
t.Xmod <- sim.comp(model,"chp5/exp-decay-inp.csv")
t <- t.Xobs[,1]; Xobs <- t.Xobs[,2]; Xmod <- t.Xmod$output[,2]
# run model evaluation
exp.eval <- model.eval(t, Xobs, Xmod, "chp5/eval-decay")

# sensitivity analysis on full factorial runs
# give nom and sd
pnom=0.02; psd =0.005; plab = c("r"); runs=5
# factorial
pval <- seq((pnom -(runs-1)*psd),(pnom+ (runs-1)*psd),2*psd)
param <- list(plab=plab, pval=pval, pnom=pnom, fact=TRUE)
# multiple run
t.X <- sim.mruns(model,"chp5/exp-sens-inp.csv", param)
# sensitivity
s.y <- sens(t.X$output, param, "chp5/exp-sens-plots-fact-out")

# sensitivity analysis on Monte Carlo runs 
pnom=0.02; psd =0.005; plab = c("r"); runs=10
pval <- round(c(pnom,rnorm(runs,pnom,psd)),5)
param <- list(plab=plab, pval=pval, pnom=pnom, fact=FALSE)
t.X <- sim.mruns(model,"chp5/exp-sens-inp.csv", param)
s.y <- sens(t.X$output, param, "chp5/exp-sens-plots-samp-out")

# gives nom lab, v in percent and levels
pnom= c(0.2,100); plab= c("r","K"); perc.v= 50; levels=3
param <- mway.fact(pnom,plab,perc.v,levels)

# use param output for simulator and msens
logis <- list(f=logistic)
t.X <- sim(logis,"chp8/logistic-100-inp csv",param,labout="sens-fact-")
s.y <- msens(t.X$output, param, "chp8/logis-100-sens-plots-factout",resp=FALSE)

# 10 runs by stratified sampling 
pnom= c(0.2,100); plab= c("r","K"); perc.v= 50; runs=10
param <- mway.samp(pnom,plab,perc.v,runs)

# use param output for simulator and msens
logis <- list(f=logistic)
t.X <- sim(logis,"chp8/logistic-100-inp csv",param,labout="sens-fact-")
s.y <- msens(t.X$output, param, "chp8/logis-100-sens-plots-factout",resp=FALSE)